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The evaluation of the number of attractors in Kauffman networks by Samuelsson and Troein 
is generalized to critical networks with one input per node and to networks with two inputs per 
node and different probability distributions for update functions. A connection is made between 
the terms occurring in the calculation and between the more graphic concepts of frozen, nonfrozcn 
and relevant nodes, and relevant components. Based on this understanding, a phenomenological 
argument is given that reproduces the dependence of the attractor numbers on system size. 
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1. INTRODUCTION 

Boolean networks are often used as generic models for 
the dynamics of complex systems of interacting entities, 
such as social and economic networks, neural networks, 
and gene or protein interaction networks [1]. The sim- 
plest and most widely studied of these models was intro- 
duced in 1969 by Kauffman [2] as a model for gene reg- 
ulation. The system consists of N nodes, each of which 
receives input from K randomly chosen other nodes. The 
network is updated synchronously, the state of a node at 
time step t being a Boolean function of the states of the K 
input nodes at the previous time step, t—1. The Boolean 
updating functions are randomly assigned to every node 
in the network, and together with the connectivity pat- 
tern they define the realization of the network. For any 
initial condition, the network eventually settles on a pe- 
riodic attractor. Thus the number and the lengths of 
the attractors are important features of the networks. 
Of special interest are critical networks, which lie at the 
boundary between a frozen phase and a chaotic phase 
[3, 4]. In the frozen phase, a perturbation at one node 
propagates during one time step on an average to less 
than one node, and the attractor lengths remain finite 
in the limit N — ► oo. In the chaotic phase, the differ- 
ence between two almost identical states increases ex- 
ponentially fast, because a perturbation propagates on 
an average to more than one node during one time step 
[5]. Based on computer simulations, the mean attractor 
number of critical K = 2 Kauffman networks with a con- 
stant probability distribution for the 16 possible updat- 
ing functions was once believed to scale as y/N [2] . With 
increasing computer power, a faster increase was seen 
(linear in [6], "faster than linear" in [7], stretched expo- 
nential in [8, 9]). Then, in a beautiful analytical study, 
Samuelsson and Troein [10] have proven that the num- 
ber of attractors grows indeed faster than any power law 
with the network size N. A proof that the number and 
length of attractors of critical K = 1 networks increases 
faster than any power law was published some time later 
[11]. These two proofs, although they apply to closely 
related systems, are conceptually different. The latter 
derives structural properties of the relevant part of the 
networks, and obtains from there a lower bound for the 



number of attractors. In contrast, in [10] the mean num- 
ber of attractors is obtained by a direct calculation that 
uses the saddle-point approximation, and which yields 
no graphic understanding of how the attractor numbers 
arise. 

It is the purpose of the present article to show how the 
attractor numbers arise in terms of the relevant parts 
of the networks. To this aim, the method chosen by 
Samuelsson and Troein is in the next section applied to 
the critical K = 1 networks, for which an intuitive un- 
derstanding already exists. The dependence of attractor 
numbers on system size N can indeed be reproduced by 
phenomenological arguments based on this understand- 
ing. In section 3, it is shown that these networks are 
similar in many respects to critical K = 2 networks, of 
which a more general class than usual will be defined. 
Applying the calculation to this more general class leads 
eventually to a phenomenological argument that repro- 
duces the dependence of attractor numbers on system 
size. 



2. CRITICAL NETWORKS WITH ONE INPUT 
PER NODE 

Let us first consider critical networks with connectivity 
K = 1. A random network with one input per node 
is critical if among the four possible Boolean functions 
only the two nonfrozcn ones, "copy" and "invert", are 
chosen. In [12] and [11], exact results for the topology 
of k = 1 networks are derived. The network consists 
of the order of ln(iV) unconnected components, each of 
which contains a loop of relevant nodes, and trees rooted 
in these loops. Relevant nodes are defined as those nodes 
whose state is not constant and that control at least one 
relevant element [9]. They determine the attractors of 
the system. The number of loops of size I is Poisson 
distributed with a mean l/l, if I is smaller than a cutoff 
size l c . The cutoff loop size scales as l c ~ y/~N [11, 12]. 

Following the calculation by Samuelsson and Troein 
[10], we calculate in the following the mean number of 
attractors of length L. More precisely, we calculate the 
mean number of cycles in state space. While an attractor 
is always a cycle in state space, the reverse is not neces- 
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sarily true, since an attractor must be accompanied by a 
shrinking state space volume. However, for the networks 
discussed in this paper, cycles are almost always attrac- 
tors, since the dynamics on the trees rooted in the loops 
is being slaved to the dynamics on the loops, and there- 
fore the initial states of the trees will be forgotten. For 
every network that contains trees, the number of initial 
states that leads to a given cycle is larger than the period 
of the cycle, and the cycles are attractors. 

Let (Cl)n denote the mean number of cycles in state 
space of length L, averaged over the ensemble of networks 
of size N. On a cycle of length L, the state of each node 
goes through a sequence of Is and Os of period L. Let 
us number the 2 i ~ 1 possible sequences of period L of 
the state of a node by the index j, ranging from to 
m — 1 = 2 i ~ 1 — 1, with sequence being the constant 
one. Following Samuelsson and Troein, we consider two 
sequences as identical if they can be transformed into 
each other by exchanging Is and Os. This simplifies the 
calculation a lot, since the sequence of the node from 
which a node with sequence j receives its input, can only 
be one sequence, which wc denote (j>(J). It is obtained 
from j by taking the first bit of j and moving it to the 
end of the sequence. Whether the Boolean function at a 
node is "copy" or "invert", has now become irrelevant, 
and all results obtained in this section apply therefore to 
critical K = 1 networks with a proportion p of "copy" 
functions and a proportion 1 — p of "invert" functions, 
for any value of p. 

If rij is the number of nodes that have the sequence j 
on a cycle of length L, and n the vector (no, . . . , n m —i), 
then 
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where (^) denotes the multinomial iV!/(no! • • • n m _i!), 
i.e., the number of different ways to assign the sequences 
to m — 1 to uq,. . . ,n m _i nodes. The factor 1/L oc- 
curs because any of the L states on the cycle could 
be the starting point, and the product is the probabil- 
ity that each node with a sequence j is connected to a 
node with the sequence 4>{j)- For sufficiently large N, 
all rij will be large, and we can apply Stirling's formula 
rij\ = (rij/e) nj ^/2nnj. Transforming the variables from 
n to x = n/N, we can replace the sum with an integral 
and obtain 



dx- , . (2) 

nm—l , 
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Itegration space is limited by the condition ^ . xj = 1 
(with all Xj > 0). The integral is evaluated using 
the saddle-point approximation, which becomes exact 
in the thermodynamic limit N — > oo. The maximum 
of the expression ^2jXjln(x^(j)/xj) is obtained when 
x 4>(j) = x j f° r &11 j- This means that the space of se- 
quences j is decomposed into permutation sets of the type 



{j, (f>(j), 4>(4>(j)), . . • }, with all members of a set occurring 
equally often at the saddle point. This can be understood 
from the topological structure of K = 1 networks. All 
nodes that are on the same component, must undergo a 
sequence belonging to the same set, while different com- 
ponents are independent from each other. Furthermore, 
on a loop or an infinitely long line of nodes, every mem- 
ber of the set occurs equally often, since between nodes 
with identical sequences, there must be nodes with all 
the other sequences from the set. The deviation from 
x <j>(j) = x j evaluated below comes from the fact that the 
branches of the trees have a finite length, which is gen- 
erally not a multiple of the set size. 

Let the index h count the permutation sets, with h = 
0, . . . , Hl — 1. Let p'l be the set with index h, which 
has 1/9^1 members. In order to perform the saddle point 
integration, we make a transformation of variables within 
each set, defining Zh = J2jep h x h an< ^ ^ x o = x j~ z hl\p\\i 
with J2jep h = 0- Only — 1 of all the Sxj within 
a set are independent. 

Expanding to second order in the Sxj, we obtain for 
the exponent in (2) 
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In the last equation, terms containing (Sxj) 2 vanish in 
the limit N — > oo, since the saddle-point integration gives 
contributions only from values Sxj of the order of 1/ \/~N. 

The integral over the Sxj can be performed by using 
the variables (Sx^u) — Sxj), leading to 




Integration space is given by the condition ^2 h Zh = 1 
(with all z h > 0). 

Let us now interpret the ./V-depcndcncc in this result. 
To this purpose, we derive the number of attractors of 
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length L from the known topological properties of K = 1 
networks. As mentioned above, the network consists of 
the order of In AT components, each of which contains a 
loop and trees rooted in the loops. The cutoff in loop size 
is l c ~ y/N. The expected number of states on a loop of 
a randomly chosen size I that belong to a cycle of length 
L is 



(5) 



the first factor being the probability that I is a multiple 
of \Pl\, and the second factor being the number of states 
in the set number h. As mentioned above, the number 
ni of loops of size I is Poisson distributed with a mean 
1/7, leading to 



L/N 
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JT e (H L -l)/l = e (H L -l)Jl" dl/l 
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e (H L -l)lnVN _ N (H L -l)/2 _ 



(6) 



The mean number of attractors of length L scales as the 
number of relevant nodes, y/N (which is proportional to 
the number of nodes in the largest loop), to the power 
Hl — 1. We will see below that an equivalent statement 
can be made for the K = 2 critical networks. In order to 
obtain also the ^-dependent prefactor in Eq. (4), the full 
probability distribution of the number of loops of a size 
of the order of l c would have to be taken into account in 
calculation (6), instead of simply integrating up to l c . 

Let us conclude this section by discussing the impli- 
cations of the fact that we do not discriminate between 
sequences that can be transformed into each other by ex- 
changing Is and Os. The numbers and the periods of the 
attractors are determined by the loops in the network. 
We call a loop "even" if it contains an even number of 
"invert" functions, and "odd" if it contains an odd num- 
ber of "invert" functions. The state of an even loop of 
size I is the same after I updates, while the state of an 
odd loop of size I is inverted after I updates. If I is a 
prime number, the period of a cycle on an odd loop is 21, 
with the second half of the cycle being obtained from the 
first half by exchanging Os and Is. However, our rules de- 
fined above (and in [10]) classify this as a cycle of period 
I, since assigning only the first half of a sequence to the 
nodes on the loop, makes a contribution to Eq. (1) if Z is 
a multiple of L. Furthermore, exchanging the Is and Os 
on a component does not lead to a new cycle according 
to our calculation, but in reality this doubles the number 
of cycles. 

Repeating calculation (6) for a system with only 
"copy" functions and by discriminating sequences that 
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TABLE I: The 16 update functions for nodes with 2 inputs. 
The first column lists the 4 possible states of the two inputs, 
the other columns represent one update function each, falling 
into four classes. 



can be transformed into each other by exchanging Is and 
0s, the result remains the same, but with Hl now count- 
ing the true number of invariant sets. 

For a system that contains also "invert" functions, the 
calculation becomes more complicated, since the mean 
number of states of a loop belonging to a cycle of length 
L is no longer given by Eq. (5). Let Hl again count the 
number of true invariant sets. The probability that a loop 
has a given cycle is now l/2|p^| if the second half of the 
cycle is not obtained from the first half by exchanging 
Is and 0s. Otherwise, the probability is 3/2|p^|. The 
mean number of states on a loop that belong to a cycle 
of length L is therefore Hl/2 for odd L and Hl/2 + Hl/2 
for even L, and these two expressions replace the Hl in 
the exponent in (6) for odd and even L respectively. 



3. A GENERAL CLASS OF CRITICAL K = 2 
NETWORKS 

Now, let us consider K = 2 networks, where each node 
has 2 randomly chosen inputs. The 16 possible update 
functions are shown in table I. 

The update functions fall into four classes [5]. In the 
first class, denoted by J- ', are the frozen functions, where 
the output is fixed irrespectively of the input. The class 
C\ contains those functions that depend only on one of 
the two inputs, but not on the other one. The class Ci 
contains the remaining canalyzing functions, where one 
state of each input fixes the output. The class 1Z contains 
the two reversible update functions, where the output is 
changed whenever one of the inputs is changed. Critical 
networks are those where a change in one node propa- 
gates to one other node on an average. A change propa- 
gates with probability 1/2 to a node that has a canalyzing 
update function C\ or C2 , with probability zero to a node 
that has a frozen update function, and with probability 
1 to a node that has a reversible update function. Con- 
sequently, if the frozen and reversible update functions 
are chosen with equal probability, the network is critical. 
Usually, only those models are considered where all 16 
update functions receive equal weight. We now consider 
the larger set of models where the frozen and reversible 
update functions are chosen with equal probability, and 
where the remaining probability is divided between the 
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C\ and C2 functions. Those networks that contain only 
C\ functions are different from the remaining ones. Since 
all nodes respond only to one input, the link to the sec- 
ond input can be cut, and we are left with a critical 
K = 1 network, which was discussed in the previous sec- 
tion. We shall see below that all the other models, where 
the weight of the C\ functions is smaller than 1, fall into 
the same class, where the number of attractors is given 
by the expression derived in [10] and reproduced below. 

For all these critical K = 2 networks, the mean number 
of attractors of length L is given by the expression [10] 
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with {Pl)\ k being the probability that a node that has 
the input sequences l\ and k generates the output se- 
quence j. This expression is the obvious generalization 
of Eq. (1) to two inputs per node. Using again Stirling's 
formula and replacing the sum with an integral, this leads 
to the generalization of Eq. (2), see [10] 
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For a network with only C\ functions, this reduces im- 
mediately to Eq. (1). The exponent has its maximum at 
zero, and this value is reached only if [10] 
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This condition is satisfied for xq = 1. For a network 
with only C\ functions, it is more generally satisfied for 
x <t>(i) — x i (f° r a U ^ 0r a ^ other critical networks, 
there exists only the maximum at xq = 1. This is shown 
as follows: Eq. (9) can be transformed into 

5>« = ^4(P £ )' 00 + 2 J] wAPl)% 
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j>0 j,k,j>0 

< x (l - x ) + (1 - x j 2 = 1 - x . (10) 

Here, we have used (Pl)qo = for j > and 
J2j>o(PL)io = 1/2 for all considered models. The in- 
equality becomes an equality only if xq = 1, or if 
Y,i,k,j>o x i x k{PL)) k = (1 - ^o) 2 - The latter condition 
is satisfied if and only if all (Pl)^ with j, k > and 
Xj , x k > vanish. They cannot vanish if there are frozen 
update functions. They do vanish if there are only C\ 
update functions. It remains to be shown that they can- 
not vanish for a system containing C2 functions, but no 



frozen functions. Assume Xi > 0. Since a node with 
two input sequences i has nonconstant output sequences 
(one of which we denote by ii) with a positive probabil- 
ity, there occurs a term XiXa(Ph)\ a- Now, the sequences 
i and ii taken together, have only 2 out of the 4 possible 
combinations of 2 bits. However, among the C2 functions 
there are functions that yield a constant output if the 
input is i and ii. Therefore even in a C2 network, not all 
(Pl)% with j, k > and Xj,x k > vanish. We thus have 
shown that all considered critical K = 2 networks satisfy 
(9) only at xq = 1. For large A, only small values Xj 
(for j > 0) contribute to the integral in (8), and a Taylor 
expansion in the Xj (for j > 0) gives [10] 



(11) 



with 
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where (A> L ) jk = (P L )j fc 



5(^X0 + <W(i))- F ° r a ^! 

network, the matrix (A l L ) vanishes, and we obtain again 
(2). The maximum of /(x) is obtained when x^us = x t 
for all i. At this maximum, the first and second term 
vanish, and the third term is of the form Ax 3 . Con- 
sequently, only values Xi (with i > 0) up to the order 
TV -1 / 3 contribute to (Cl)n- This means that the pro- 
portion of nodes that are not frozen on an attractor is 
of the order TV -1 / 3 , and the total number of nonfrozen 
nodes is of the order A 2 / 3 . This is in contrast to the 
critical K = 1 network, where a nonvanishing proportion 
of nodes is nonfrozen. Changing the variables again to 
z h = Yjiapl Xu an< ^ = X% ~ Zh /\Pi\' tne integration 
over the Sxi gives now 
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with {B h L ) jk = (P' L ) h ]k - \{5, ]h + 5 kh ), and with (P' L )% 
being the probabity that the output sequence belongs to 
set h if the input sequences belong to the sets j and k. 
Introducing a new variable = z^A 1 / 3 , we obtain an 
additional factor Af-t^- 1 )/ 6 , and the mean number of 
cycles of length L becomes [10] 



Yl h>o 



ir \ ~ 1 N " L3 1 TT J dyh 

(14) 

While integration space for the z% was restricted by the 
condition Zh = 1 — xq, there is no constraint for the 
Vh- 
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With the understanding gained from the K = 1 criti- 
cal networks, we can interpret the calculation as follows. 
The difference between K = 1 and K = 2 critical net- 
works comes from the fact that in the K = 2 networks 
only the fraction iV -1 / 3 of nodes is nonfrozen. This mod- 
ifies the exponent of N in the final result, and this leads 
to the different form of the Zh integration. Both types 
of networks have in common that the main contribution 
to the integral comes from the neighborhood of the sub- 
space satifsying Xm{\ = Xi for all i. This means that the 
majority of nonfrozen nodes receive input from one non- 
frozen node, the other input being frozen. The nonfrozen 
part of a K = 2 critical network resembles therefore a 
K = 1 critical network. The proportion of nonfrozen 
nodes receiving input from two nonfrozen nodes, cannot 
be larger than of the order iV -1 / 3 , since the Sxi are of 
the order A^ -1 / 3 . Thus, the nonfrozen part of a K = 2 
critical network differs from a K = 1 critical network by 
a proportion N -1 ' 3 of nonfrozen nodes having two non- 
frozen inputs. Apparently, this difference does not affect 
the scaling of {Cl)n with N, but only the L-dependent 
prefactor. If the number of relevant nodes scales as N 1 ^ 3 
(as is numerically found in [7]), the law 

means that the mean number of attractors of length L 
scales as the number of relevant nodes, N 1 / 3 (which is 
proportional to the number of nodes in the largest com- 
ponent), to the power Hl — 1. It can be obtained by a 
phenomenological argument similar to the one used in the 
previous section. Assume there are N 1 / 3 relevant nodes 
arranged in ~ In AT components, with the number of com- 
ponents of size I being Poisson distributed with a mean 
l/l. Since at most the proportion N^ 1 /' 3 of nonfrozen 
nodes have two nonfrozen inputs, only a finite number 
of relevant nodes have two nonfrozen inputs, and all rel- 
evant components apart from a finite number are cycles 
without additional links, just as for the K = 1 critical 
network. The mean number of states on a component of 
size I that belong to a cycle of length L is therefore in 
the limit of large N 

just as for the K = 1 critical network. If the number n; 
of relevant components of size I < l c ~ N 1 / 3 is Poisson 
distributed with a mean l/l, we obtain 

<^>* - e n ( ^F-^') 

{ni}l<la \ / 
_ ^ffi-lJlnJV 1 ' 3 = N (H L -l)/3 _ / 15 ) 



4. CONCLUSIONS 



In this paper, we have considered the mean number 
of attractors of length L for critical K = 1 and K = 2 
networks. We have applied the method by Samuelsson 
and Troein [10] and have interpreted the results in terms 
of the topological properties of the nonfrozen part of the 
network. For the K = 1 networks, the dependence of 
the number of attractors of length L on the system size 
N, (Cl)n ~ N i ^ Hl ^ 1 ^ 2 could be understood as resulting 
from the network containing of the order of N 1 / 2 relevant 
nodes arranged in ~ In A" components, with the number 
of components of size I being Poisson distributed with a 
mean l/l. The nonrelevant nodes sit in trees rooted in 
the loops. 

Then, we could show that all K = 2 critical net- 
works can be treated by the same calculation. Only 
for networks consisting only of C\ functions, the step 
from Eq. (13) to Eq. (14) cannot be made, since the 
matrix (Bfyjk vanishes in this case. Ci-networks are 
in fact K = 1 critical networks, and Eq. (13) is iden- 
tical to Eq. (4) in this case. All the other K = 2 net- 
works show the same dependence of attractor numbers 
on system size, with only the L-dcpendent prefactor be- 
ing different (because the matrix {B^jk is different for 
a different choice of weights for the update functions). 
We saw that only the proportion N~ x / 3 of nodes is non- 
frozen, and that almost all nonfrozen nodes depend only 
on one nonfrozen input. The nonfrozen part of critical 
K = 2 networks resembles strongly a K = 1 critical net- 
work, and by analogy to the K = 1 critical network we 
concluded that the scaling with A^ of the number of at- 
tractors of length L in K = 2 critical Boolean networks 
can be understood as resulting from the network being 
composed of the order of N 1 / 3 relevant nodes arranged 
in ~ In AT components, with the number of components 
of size I being Poisson distributed with a mean l/l. Only 
the proportion N" 1 ^ 3 of all nodes is not frozen, and those 
nonfrozen nodes that are not relevant sit in trees rooted 
in relevant components. 

The calculation indicates that there are several quan- 
tities that show a scaling with N. Among these are 
the number of nonfrozen nodes, the number of relevant 
nodes, and the number of nonfrozen nodes with two 
nonfrozen inputs. Evaluating these scaling properties 
to more detail could be the next step in understanding 
Kauffman networks. 

I thank F. Greil and T. Mihaljev for comments on the 
manuscript. 
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